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Multiple-type branching processes that model the spread of infectious diseases are investigated. In 
these stochastic processes, the disease goes through multiple stages before it eventually disappears. 
We mostly focus on the critical multistage Susceptible-Infected-Recovered (SIR) infection process. In 
the infinite population limit, we compute the outbreak size distributions and show that asymptotic 
results apply to more general multiple-type critical branching processes. Finally using heuristic 
arguments and simulations we establish scaling laws for a multistage SIR model in a finite population. 
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I. INTRODUCTION 

A critical branching process (CBP) is the fundamen- 
tal stochastic process that has innumerable apphcations 
[IH5]. Although the adjective "critical" suggests that 
the CBP is a peculiar branching process which should 
require the tuning of parameters, the CBP is arguably 
more important than its generic sub-critical and super- 
critical brethren. Indeed, the sub-critical branching pro- 
cess quickly leads to extinction while the super-critical 
branching process may result in an unlimited growth, 
so whenever we are seeking a framework for phenomena 
which are driven by branching and are maintained in a 
stationary state, the CBP provides the proper setting. 

The simplest branching process involves duplication 
and death, and hence called a birth-death process. Sym- 
bolically 

A + A 



A 



For the critical birth-death process the probabilities of 
these two events are equal. One can think on As be- 
ing cells which either divide or die. Another interpreta- 
tion is to identify As with infected individuals who ei- 
ther die or transmit infection to other individuals. The 
branching process then becomes a stochastic version of 
the Susceptible-Infected-Recovered (SIR) infection pro- 
cess. According to this model [MS], the population con- 
sists of susceptible, infected, and recovered individuals, 
and the infection spreads by contact between infected and 
susceptible members of the community, while infected in- 
dividuals recover (and become immune), or equivalently 
they are removed from the system. The branching pro- 
cess also provides the simplest model of biodiversity [10] 
(see for a survey of biological literature on macroevo- 
lution modeling), with duplication and death processes 
corresponding to speciations and extinctions. The CBP 
is neutral as it assumes that speciations and extinctions 
occur with equal probability. In the context of infectious 
diseases the CBPs naturally arise as a balance between 
human efforts leading to the reduction of the infection 



rate [7] and the natural evolution that increases the in- 
fection rate of diseases hovering below the threshold [12] . 

In numerous applications, it is desirable to general- 
ize the CBP to include more than one type of elemental 
building blocks. In cell biology, we may have two pop- 
ulations of cells, progenitor A cells that can divide and 
differentiate into B cells. These two cell populations can 
be (on average) in a stationary state (homeostasis), if 
they evolve according to the two-type CBP. The simplest 
two-type CBP is represented by the scheme 



A + A 
A 



A- 



B + B 
A 



B 



where all steps are equiprobable. 

Branching processes can be studied in discrete time 
when all infected individuals divide or die synchronously, 
or in continuous time when each infected individual inde- 
pendently divides or dies. For a critical process the birth 
and death probabilities or rates are equal. Continuous- 
time two-type birth-death branching processes are solv- 
able, namely the time-dependent sizes of populations of 
A and B individuals have been analytically determined 
Specializing the results of Refs. [T3] to the above 
simplest two-type critical CBP leads to simplifications, 
yet the results are still cumbersome. In this paper we 
are concerned with different questions, namely we are in- 
terested in time-independent quantities summarizing the 
outbreak size distributions. 

The rest of this paper is organized as follows. In Sec. [IT] 
we consider the simplest multiple-type critical branch- 
ing process and derive, in the infinite population limit, 
the outbreak size distributions. In Sec. IIIII we show the 
universality of the asymptotic results for the outbreak 
size distributions, namely that at the epidemic threshold 
these distributions have algebraic tails which are deter- 
mined only by the number of types (but not the detail 
of the model). In Sec. IV we analyze finite populations 
and obtain finite-size scaling laws for the sizes of the out- 
breaks. In Sec. [Vjwe give a short discussion. Finally, in 
an appendix, some of the details of the analysis of a gen- 
erating function are presented. 
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II. CRITICAL BIRTH-DEATH PROCESS 

In this section we consider the infinite population hmit. 
We shall always assume that the system begins with a sin- 
gle A, representing an infected individual. In the multi- 
type case, an A can turn into B, a B into C, and so 
on. The probability for each individual to duplicate or 
disappear is the same, which makes our process critical. 
Consequently, all individuals eventually disappear with 
probability one and the epidemic stops. We want to de- 
termine the outbreak size distribution. 



A. Single Type Critical Birth-Death Process 

Let An be the probability that exactly n individuals 
catch the infection before the epidemic is over. We have 
Ai — ^. Further, A2 — \A\ since at the first step a 
new individual must get infected and then both must die 
without spreading infection. Proceeding along these lines 
we arrive at the recurrence 



1 



(1) 



reflecting that the first infection event results in two in- 
dependent infection processes [2]. A solution to ([I]) is 
found by introducing the generating function 



n>l 



(2) 



which converts ^ into a quadratic equation 2A — A'^ + z 
whose solution reads 



A{z) = 1 - Vl- z 
Expanding A{z) in powers of z we find 



(3) 



1 T{n-l] 



1 



-3/2 



r(n + i) ^ 

In particular, the probabilities An are given by 
1 1 1 5 7 21 33 429 715 2431 



(4) 



2' 8' 16' 128' 256' 1024' 2048' 32768' 65536' 262144 

for n = 1, . . . , 10. Note that the above process can also 
be interpreted as a discrete time symmetric random walk 
on the integers, where An corresponds to the probability 
that the first return to the origin occurs at time 2n [T] . 



B. Two- Type Critical Birth-Death Process 

In this case the size distribution of the outbreaks of A 
types is the same as for the single type CBP. Consider 



now the outbreaks of B types. Denote by B„ the proba- 
bility that exactly n individuals of type B are born before 
the outbreak is over. These probabilities satisfy 



BiBi 



^ An 

2 



(5) 



To understand this recurrence it suffices to consider what 
happens to the initial A individual. One possibility is 
that the initial A turns into a B (this happens with prob- 
ability 1/2), and then this B produces n individuals of 
type B with probability An- Combining we get the sec- 
ond term |j4„ on the right-hand side of Eq. Oth- 
erwise the initial A produces two As (this happens with 
probability 1/2), and then the probability that these two 
As will produce exactly n individuals of type B is just 
the convolution J2i+j=7i BiBj. Overall this gives the first 
term on the right-hand side of Eq. ^ . 

To tackle Eq. ^ we again proceed by introducing the 
generating function 



S(z) = Y.Bn-^ 



(6) 



which converts ^ into a quadratic equation 223 — 'B^+A 
whose solution reads 

a3(z) = l-(l-z)V4 (7) 
Expanding 23(2:) in powers of z we find 

1 r(n-i) 1 



Bn 



4r(|) r(n + i) 4r(f) 



-5/4 



(8) 



C. Multiple- Type Critical Birth-Death Process 



Generally for multiple-type CBP we seemingly have 



2 ^ A^^Af^ + -^AlP-^^ 



(9) 



The corresponding generating functions A'^p and A'^p'^ 
are related via 



2^1 = [A'-P^f+A^P-^^ 
We already know that 

yiw = yi = 1 - (1 - z)i/2 

= S = l-(l-z)l/4 



(10) 



A general solution to ( 10 1 is 



from which 

AiP^ 



^(f) = 1 - (1 - zy/^' 

r{n~2-P) 



(11) 



1 



2pr(i-2-p) r(n + i) 2pr(i-2-p) 

. ^ . . (^^^ 

One can verify (111 by induction. It also follows from the 



more general result of the following subsection II D 
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D. Joint Outbreak Size Distribution 

The joint outbreak size distribution first becomes rel- 
evant for the two-type CBP. In this setting we want to 
compute the probability Pm,n that exactly m individuals 
of type A and n individuals of type B are born before 
the outbreak is over. 

The disappearance of an A results in the birth of a _B, 
and therefore Pm.n = whenever m > n. To compute 
Pm,n = in the interesting m < n range, we shall use 
the important abelian property of the model: The com- 
position of the system at the moments of birth of new Bs 
is irrelevant, so one can assume that the process starts 
with m individuals of type B (and no A individuals). We 
have Pm,m = Am{Ai)™ where the first term assures that 
exactly m individuals of type A ever exist and the second 
term gives the probability that all m individuals of type 
B eventually disappear without duplication. Further we 
have Pm.m+i ~ Amm{Ai)™-~^ A2, and generally 



Pjn^n — ^ ^ A{-^ . . . A^ 

Using the generating function 

n>m>l 



(13) 



(14) 



and ( 13 1 we get 

y{x,y)=Y,A^x^[A{y)Y' 

m>l 

= A{xA{y)) 
= 1 - - =^My) 

= 1 - \Jl-x + - y 

Specializing a; or y to unity wc recover previous results 
^ and 0: 

y{x,y^l)=A{x), 7(2; = l,y) -S(y) (15) 

This essentially provides another derivation of Eq. ([7|. 

From the generating function ( 14 ) we can express the 
probabilities more explicitly as 



P„ 



(2m -3)!! 
' 2"'+"n! '^"'"-'"+1 '-"-"^ 

where the integer Sn,k for n > k > 1 are given by the 
recursion 



Sn,k = Sn-l,k + {n + k- 3)s„_i^fe_ 



(17) 



with Si^i = 1; otherwise we consider Sn,k to be zero. The 
first few values of s„ are as follows 



/ 1 \ 

1 1 

1 3 3 

1 6 15 15 

1 10 45 105 105 

V 1 15 105 420 945 945 / 



(18) 



For example, the above formulas lead to the following 
probabilities Pm.n for m,n < 6 



/ 1 J_ J_ 

' 4 16 3,2 

I2 









V 



5 



21 

2048 



7 
5J,2 

1(^4 4(^6 



128 5J2 2(^48 2y^ 



8 



2048 






2C;^8 

8192 32768 


" 65536 



16^4 



(19) 



One can recognize the first few rows of the matrix Pm.n- 
For the first three rows 



Pun = 
P2,n = 
Pz.n = 



(2n-2)! 1 
(n - l)!n! 2^^ 
(2n-2)! 1 



{n-l)\n\ 22"+i 
(2n-4)! 3 
(n-3)!n! 22"+i 



n>2 



n > 3 



(20) 



These results are straightforwardly verified by induction. 
One can also identify the diagonal elements 



Pn. 



1 



(2n-2)! 
(n- l)!n! 23«-i 



(21) 



We haven't succeeded, however, in finding an explicit ex- 
pression for Pm.n in the general case. 

Consider now the three-type CBP. We again limit our- 
selves with the simplest three-type CBP which is repre- 
sented by the scheme 



A + A 



A- 



B + B 



B 



C + C 



■ c- 



where all events occur with probabilities |. The prob- 
ability Pm.n,e that exactly (m, n, £) individuals of types 
{A, B, C) are born before the outbreak is over is given by 
an obvious generalization of Eq. : 



Pm,n,t — Am ^ ^ • • ■ ^'i,„ ^ ^ ^ji ■ • ■ Aj^ 

Jl + -+j„=e 



ii + . . . + im—n 



Converting this recurrence into an equation for the gen- 
erating function 



y{x,y,z)^ Prr 



l>n>m>l 



we find 

y{x,y,z) = A{xA{yA{z))) 



= I — yl — x + xyl — y + yy/l — z 

Generally for the type branching process the gener- 
ating function 



?(zi,...,Zp) = Y 



nil , .. . .rrip-^i 



mi rrip 

-C-i ... 4-p 



mp>...>mi>l 
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which encodes the joint distribution Pmi....,m can be ex- 
pressed through the generating function A of the single 
type process as 

^(zi, . . . , zp) = A{ziA{z2A{. . . A{zp_iAizp)) ...))) 

(22) 

SpeciaHzing y{zi, . . . , Zp) to zi = . . . = Zp_i = 1, Zp — z 
we must reproduce the generating function A'^p^z). We 
actually get an identity 



A^P\z) ^A{A{...A{A{z))...)) (23) 
Using ( pS] ) in conjunction with ([3| we indeed recover (111. 



III. MORE GENERAL PROCESSES 

Let us first generalize the results of the previous section 
to non-critical brith-death processes. For a birth-death 
branching process we set the recovery rate to unity and 
denote by a the infection rate. The process is called sub- 
critical for a < 1, critical for a = 1, and supercritical 
for a > 1. While for a < 1 the outbreak is finite with 
probability one; for a supercritical process (a > 1) the 
outbreak never stops with probability 1 — 1/a. We can 
still calculate the outbreak size distribution by condition- 
ing on finite outbreaks. Let A„ be the probability that 
a single individual results in a finite outbreak where ex- 
actly n individuals are infected during the epidemic. We 
can write the following recursion 

1 + a 1 + a 

i+j—n 

which, for the generating function takes the form 

{l + a)A^aA^ + z (25) 
and its solution is given by 



2a 



The probability of a finite outbreak is of course 
■^(^ = 1) = { l/a foi a > 1 



(27) 



and the average outbreak size in case of a finite outbreak 



(n)=5]nA„=yi'(z = l) = ^4^ (28) 

By expanding A in powers of z, we obtain the proba- 
bility of a finite outbreak of size n 



1 



4a 



(1 



"r(n-j) 
r(n + 1) 



(29) 



Asymptotically, for large n we have a power law decay 
1 -I- a 



Arn 



with an exponential cutoff around 

^= [21og(l + a)-log4a]^^ 



(30) 



(31) 



The cutoff diverges as a — > 1, and for the critical process 
a = 1 we find a pure power law decay Q. 

More generally, consider a branching process where 
each individual is replaced by k individuals at rate Uk. 
Equivalently, we can consider a discrete time branching 
process, where at every time step each individual is in- 
dependently replaced by k individuals with probability 
Pk = OikI ^oin- In this case the governing equation for 
the outbreak distribution (which is essentially the back- 
ward Kolmogorov equation) is a simple generalization of 
Eq. (25), namely 



(32) 



fc>0 



k>l 



An explicit solution becomes less accessible for larger 
number of possible offspring, but for example the crit- 
ical process with three offspring (ao = 2/3, — 1/3; 
other rates are zero) with a single initial individual is 
still tractable 



yi = 2 sin ■ 



(33) 



Note that in this case it is impossible to find a cluster 
with even number of individuals: = 0. The first few 
nonzero values of A„ are 



, 2 ^ 8 . 32 ^ 

^^ = 3'^^ -81'^^- 729'^^ 



512 , 28160 



19683' 



1594323 



From the behavior of yi(z) around z = 1 we find that for 
large odd n the asymptotic is 



1 



-3/2 



(34) 



The exponent is again equal to 3/2. Using (32) one can 
show that this exponent is universal for any (single-type) 
branching process. We state this for the equivalent dis- 
crete time version where the probability of a individual 
having k offspring is pk = oik/ 1]„ for A: 0, 1, 2, ... : 

Consider a discrete time branching process with a single 
initial individual and an offspring probability function pk ■ 
Suppose that the process is critical, "^kpk = I, and has 
finite variance = ^k^Pk — 1 < oo. Then asymptoti- 
cally 1 — A ^ ■y/c(l — z) for z — > 1, and 



for n — >■ oo, with c — 2po/cr'^ 



-3/2 



(35) 
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To prove this assertion we recall that A{z) is analytic for 
\z\ < 1 and A{1) = 1. We write A{z) = 1 - 'B{z), where 
lini^_j.i 2?(z) = 0, and plug this expression into (32). Us- 
ing the conditions for the zeroth and first moments of pk , 
we find 



1-03 = Y.p,{l-'B)'-po{l~z) 

+ OCB') ~ Po{l - z) 



k>0 



= 1-03 + — '--3. 



which gives the leading asymptotic behavior 

S = v/c(l - z) + 0(1 - z) (36) 



with c ~ 2po ja^ . Expanding \/\ — z, or equivalently us- 
ing a Tauberian theorem [l] , we arrive at the asymptotic 



(35) for Am which completes the proof. 



The asymptotic behavior of multi-type processes now 
follows from recursion (|23l): 



1 - A^p\z) = 1 - yi(yi(f-i)(z)) - c^i'^[i - yi(p-i)(z)]i/2 

^c^'^c^'^[l-A^P-^\z)Y/^ ^ ... 
^c'-^-'il-zf-' 

If at all stages there are different critical branching pro- 
cesses, then simply 



\-A(^\z)^c\'\\'^---cl'^''{l-zf 



(37) 



where all Cj are defined analogously. This corresponds to 
the asymptotic behavior for the probabilities 



An OH n 



-1-2-" 



(38) 



IV. FINITE SIZE EFFECTS 



Epidemic outbreaks obviously involve finite popula- 
tions. In other applications of branching processes the 
finiteness also plays a crucial role. Finite size effects are 
especially pronounced for the critical branching processes 
[TH - [ro| . Therefore it is important to understand how the 
basic characteristics such as the average outbreak size 
[151 116| vary with population size N. 

In a standard SIR model individuals of a population of 
size N are either susceptible to the infection (S), infected 
(I), or recovered (R), who cannot get infected again. In a 
multistage model infected individuals can be at different 
stages of infection, and for all (except the last) stages 
recovery means passing to the next stage. The recov- 
ery rates are equal to one as in the infinite-population 
limit. The infection rate, however, is only initially equal 
to unity (as the process is critical) , and it decreases with 
time since only susceptible individuals can get infected, 
and their number decreases. More precisely, the infec- 
tion rate at any time is the ratio of the total number of 
susceptible individuals to the total population size. 



The average size (n) = J2n>i epidemic out- 

break diverges in an infinite system due to the power-law 
tail ^ of the probability distribution v4„. Setting the 
upper bound at the population size N one estimates the 
average size of the epidemic outbreak 



N 



n=l 



-1/2 



N 



(39) 



This heuristic argument is actually based on the tacit 
assumption that a finite fraction of the population [as 
reflected by the upper limit in the sum in (39)] may get 
infected. It turns out that the actual average outbreak 
size is much smaller because the epidemic outbreak weak- 
ens as more individuals become infected. Let M be the 
maximal outbreak size. The same estimate as in (39) 
now gives (n) oc y/M. On the other side, the effective in- 
fection rate is equal to (A^ - M)/iV = 1 - M/N, and the 
average size of an outbreak in such sub-critical branching 
process is N/M, see (28). These two estimates should 
match, y/M oc N/M, from which M a N'^/^ thereby 



amending the naive estimate ( 39 ) to a self-consistent es- 
timate [16] 



(n) cx v/]V 



(40) 



The scaling law ( 40 ) gives the correct asymptotic behav- 



ior. This interesting result was established by Martin-L6f 
[l5] ; several more recent studies [I6H20] confirm ( 40 1 and 
analyze other finite-size effects. 

Consider now the two-type CBP in the context of epi- 
demic processes. We can think about As as individuals 
carrying light infection which can be transmitted to sus- 
ceptible individuals and evolve into deadly infection. The 
latter can be also be transmitted to susceptible individ- 
uals or it can cause death. We need to find the scaling 
behavior of the maximal outbreak sizes Ma and Mb- 
Having established AI^ and Mg, one can use the power- 
law tails ([4]) and ([8| to deduce the average outbreak sizes 



(ng) cx M' 



3/4 



(41) 



We again obtain an estimate ^/M^ (x N/Mg. A naive 
conjecture is that the maximal outbreak sizes are com- 
parable. In this situation 



Ma oc Mb oc N^'^ 

(ua) oc Vn , (ub) oc VN 



(42) 



We know, however, that different types follow different 
scaling laws, specifically A^ oc n~^^^ and i?„ oc n^^/'', 
and hence it seems much more plausible that Ma <C Mb . 
In this case we use 



N 
Mb 



N 



'Ma — oc M, 



3/4 



(43) 



The second relation in ( 43 ) is understood by noting that 
Eq. ^ with a = 1-A'I^N gives N /Mb for the average 
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FIG. 1: Mean outbreak sizes arising in simulations for a two- 
type GBP in finite populations from an average over 10^ runs. 
The slopes of the lines are given by ( 45 1 . 




FIG. 2: Mean outbreak sizes arising in simulations for a three- 
type GBP in finite populations from an average over 10^ runs. 
The slopes of the lines are given by ( 47 1 . 



size of the outbreak caused by a single B, so in the case 
when the number of the seed individuals of type B is 
equal to V-^A, we get \/WaN/Mb for the average size 
of the outbreak. From ( 43 ) we obtain the scaling of the 
maximal outbreak sizes 



Ma oc TVtt ^ Mb oc N~i 



(44) 



which is combined with (41 1 to yield the average sizes 
(tm) oc , (hb) cx (45) 



For the three- type CBP the analog of ( 43 1 reads 
N 



cx M. 

Mr 



1/2 



M 



1/2 



N 



M, 



cx M' 



3/4 



^ Mn ^ 



from which the maximal outbreak sizes scale as 

14 56 24 , , 

Ma cx iVsT , Mb cx , Mc cx N^i (46) 
and the average outbreak sizes are given by 

(ua) cx iV^r , (tib) cx , (nc) cx iV^r (47) 



The predictions ( 45 1 and (|47[) for the average outbreak 



sizes in the two- and three-type CBPs are in good agree- 
ment with simulations, see Figs, [l] and [2] 

Generally for the p— type CBP one finds that the max- 
imal and average outbreak sizes scale as 

Mq CX iV^p9/(l-2-') ^ ^^^^ ^ J^X.q (48) 

where q — 1,2, . . . ,p labels the types of individuals and 



Xr, 



l+p-i(l-2-P) 



The outbreaks of individuals of type p (those which die 
rather than differentiate) are the largest: 

MpCxiVi/[i+P"'(i-2"")l 



V. DISCUSSION 

We investigated the multistage SIR infection process 
at the epidemic threshold. In an infinitely large popu- 
lation, we computed the outbreak size distributions and 
showed that these distributions have universal algebraic 
tails. Combining these infinite-population results and 
heuristic arguments we established scaling laws for finite 
populations. Specifically, we showed that the maximal 
and average outbreak sizes scale as powers of the popula- 
tion size. For instance, for the three-type critical branch- 
ing process six non-trivial exponents, see Eqs. (46)-(47l, 



describe the maximal and average outbreak sizes of the 
three types. We numerically verified the scaling laws and 
found an excellent agreement. 

A number of extensions are worth pursuing, especially 
those which would increase our understanding of the fi- 
nite size scaling. Many of these extensions are easy to for- 
mulate, yet very challenging to analyze. For instance, one 
would like to prove the validity of the chief scaling results 
( 48 1 . Even more demanding would be to determine the 



probabilities A„(7V), Bn{N), etc. For the two- type CBP, 
for instance, the ratios of the probabilities An{N), Bn{N) 
to their infinite-population values An{oo), Bn{oo) given 
by Eqs. Q, are expected to exhibit scaling behaviors 



An{oo) 



F 



(ats/h) 



Bn{oo) 



G 



ATS/ii 



(49) 



The scaled sizes in Eq. ( 49 1 are fixed by the scaling laws 



(44 1 for the maximal outbreak sizes. 



7 



It would be also interesting to determine the duration, 
both average and maximal, of the outbreaks; for the clas- 
sical critical SIR process, this problem has been investi- 
gated in [HITnillO]- 

Finally, we stress that in this paper we limited our- 
selves to stochastic processes without any spatial struc- 
ture. Needless to say, in applications the spatial or net- 
work structure of the infected domain play an important 
role [?IVt27j . The generalizations to such settings form 
an interesting avenue for future research. 
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Appendix A: Derivation of two-type probabilities 



We have obtained the explicit expression ( 16 ) by dif- 
ferentiating the generating function 



p — 



1 



7 ' 77 . ' 



(Al) 



Here we provide some details. First, we differentiate the 
generating function (|14[) with respect to x to obtain 



(2to-3)!! (l-Vr^)" 



(A2) 



We can set a: = on the right-hand side of (|A2j) since it 
won't affect Pm,n defined by (All. Hence we should dif- 
ferentiate 9™J'U=o = ^^^^#^(l-yr^)™ with respect 
to y only. At each time we differentiate (1 — ^/l — y)™, 
the number of terms grow by one until n = rii, at which 
point the m terms are proportional to (1 — y/1 — yY with 
I = 0, . . . , m — 1. Differentiating further won't increase 
the number of terms. Therefore 



^ 2»(m-n-l-h/c)!(l-w)("+*:-i)/2 

This can be proved by induction and by using the recur- 
sion formula ( |17[ ) for Sn,k- Setting y = in the sum on 
the right-hand side, we are left with only one term (the 
smallest k term). Accounting for all factors in (Al) we 
arrive at the announced result (16). 
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